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Abstract 


Ensemble and reduced-rank approaches to prediction and assimilation rely on low-dimensional 
approximations of the estimation error covariances. Here stability properties of the fore- 
cast/analysis cycle for linear, time-independent systems are used to identify factors that cause 
the steady-state analysis error covariance to admit a low-dimensional representation. A useful 
measure of forecast/analysis cycle stability is the bound matrix, a function of the dynamics, 
observation operator and assimilation method. Upper and lower estimates for the steady-state 
analysis error covariance matrix eigenvalues are derived from the bound matrix. The estimates 
generalize to time-dependent systems. If much of the steady-state analysis error variance is 
due to a few dominant modes, the leading eigenvectors of the bound matrix approximate those 
of the steady-state analysis error covariance matrix. 

The analytical results are illustrated in two numerical examples where the Kalman filter is 
carried to steady state. The first example uses the dynamics of a generalized advection equation 
exhibiting nonmodal transient growth. Failure to observe growing modes leads to increased 
steady-state analysis error variances. Leading eigenvectors of the steady-state analysis error 
covariance matrix are well approximated by leading eigenvectors of the bound matrix. The 
second example uses the dynamics of a damped baroclinic wave model. The leading eigenvec- 
tors of a lowest-order approximation of the bound matrix are shown to approximate well the 
leading eigenvectors of the steady-state analysis error covariance matrix. 
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1. Introduction 


Accurate and realistic representations of forecast and analysis errors are key to the performance 
of atmospheric data assimilation and ensemble prediction systems. Unfortunately, forecast and 
analysis error statistics such as the covariance are impractical to compute directly for operational- 
scale problems, and approximate methods must be used. Primary sources of difficulty are the large 
number of degrees of freedom present in numerical forecast models and poor understanding of 
forecast model error, making it infeasible to store and evolve the complete error covariances. 

Data assimilation systems combine observations with a forecast model first guess to produce an 
analysis or estimate of the state of the atmosphere. If the error statistics of the forecast and of the 
observations are known, then the relative weights given to the first guess and to the observations 
can be chosen so that the analysis error is minimized. Data assimilation systems are sensitive to the 
specification of the forecast error statistics, and misspecificatiori of them increases analysis error 
levels. In current operation practice, the forecast error covariance is generally approximated by 
an analytical state-independent model depending on a number of estimated parameters (Courtier 
et al., 1998; Rabier et al., 1998; Dee and da Silva, 1999; Dee et al., 1999). Though the details 
of these error models vary, their simplicity means that realistic features, such as anisotropy and 
flow-dependent structure, are not included. A general way of including some of these features is 
through reduced-rank methods, where evolution of the error covariances or the error covariances 
themselves are replaced by low-dimensional representations (Evensen, 1994; Cohn and Todling, 
1996; Cane et al., 1996; Verlaan and Heemink, 1997). 

Ensemble prediction schemes produce an ensemble of forecasts, each starting from slightly 
different initial conditions. The quality of the ensemble forecast, particularly for short to medium- 
range forecasting, is sensitive to the structure of the perturbations used to form the ensemble of 
initial conditions. Ideally, the ensemble of initial conditions should reflect the errors present in the 
analysis. Operational ensemble forecast systems use state-dependent initial perturbations that are 
related to either past or future growing modes of the dynamics (Toth and Kalnay, 1993; Molteni 
et al., 1996). The existence of a distinction between past and future growing modes stems from 
the time-dependence and nonnormality of the tangent linear dynamics which is linearized about 
a nonlinear, time- varying trajectory. If the dynamics were time-independent and normal (a lin- 
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ear operator is normal when it has a complete set of orthogonal eigenvectors), both past and future 
growing modes would coincide with the eigenvectors or normal modes. Time-independent dynam- 
ics. for instance, the result of linearizing about a mean state or of statistical modeling, will typically 
be nonnormal and present behavior that is primarily nonmodal (Blumenthal, 1991; Whitaker and 
Sardeshmukh, 1998). The rationale for using past growing modes (approximate Lyapunov vectors) 
as initial perturbations is that they represent the analysis error due to errors in the first guess. How- 
ever, there is no guarantee that past growing modes efficiently sample the space of rapidly growing 
analysis errors. On the other hand, future growing modes (singular vectors of the dynamics), while 
presenting optimal growth, intrinsically have little to do with the analysis error. If information 
about the analysis errors is available, the singular vectors of the dynamics can be used to identify 
the analysis errors that will contribute most to forecast error (Barkmeijer et al., 1998; Swanson 
et al., 1998). In any case, since an ensemble covariance is by construction low-rank, it is useful to 
identify factors that lead to the actual error covariances having good low-dimensional approxima- 
tions. 

In this article the analysis error covariance structure is examined in the context of linear, time- 
independent dynamics and observation operators. The aim is to understand under what conditions 
the analysis error covariance admits a low-dimensional representation. Such information would 
be useful in both ensemble forecast systems and data assimilation systems, though our results in 
this idealized context can only give qualitative information. Our approach is to identify errors 
which are not efficiently reduced by the forecast/analysis cycle. If the system is observable and the 
Kalman filter, or more generally a stable data assimilation method, is used, then all analysis errors 
are eventually reduced, that is to say, the forecast/analysis cycle is exponentially stable. A linear 
system and its measurements is completely observable when it is possible to reconstruct uniquely 
the state of the system from past measurements (Cohn and Dee, 1988). However, it may be the 
case that the damping of analysis errors by the forecast/analysis cycle is slow or that analysis errors 
may be amplified nonmodally on transient time-scales before being reduced. Only in special cases 
is the behavior of forecast/analysis cycle purely modal, for instance if the dynamics is normal, all 
variables are observed, and the observational and forecast model errors are homogeneous (Daley 
and Menard, 1993). 

Our motivation is complementary to that of Swanson et al. (1998), who studied the analysis er- 
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ror structure of four-dimensional variational assimilation applied to a nonlinear quasi-geostrophic 
model. There, in the perfect model setting with all state variables observed, it was found that in 
the limit of a long assimilation period, analysis errors projected primarily onto Lyapunov vectors. 
These Lyapunov vectors were found to project weakly onto the leading right singular vectors of 
the dynamics, although this part of the analysis error was responsible for most of the forecast error. 
Distinguishing features of the present work are treatment of the interaction of observing pattern and 
dynamics, and analytical results that include model error but do not require its complete details. 

We make extensive use of the fact that the steady-state analysis error covariance matrix sat- 
isfies a discrete algebraic Lyapunov equation (DALE) in the case of time-independent dynamics, 
observation operators, forecast model error and observation error covariances. Investigation of 
the DALE permits the identification and estimation of the dominant part (leading eigenvalues and 
eigenvectors) of the steady-state analysis error covariance and understanding of its dependence on 
the dynamics and observing pattern. The transient and asymptotic behavior of the analysis errors, 
as measured by the bound matrix, provides rigorous estimates for the eigenvalues and eigenvectors 
of the steady-state analysis error covariance. A bound matrix giving time-dependent estimates for 
the time-dependent analysis error covariance can also be defined when the dynamics, observing 
pattern and error covariance sources vary in time. 

The organization of this paper is as follows. In Section 2 the generic linear, time-independent 
forecast and analysis system is described. In Section 3 the bound matrix is defined and employed 
to investigate the steady-state analysis error covariance matrix. Sections 4 and 5 illustrate the 
results of Section 3 by applying the Kalman filter to two simple dynamical models: a generalized 
advection equation and a model for damped baroclinic waves. Finally, a summary and conclusions 
are given in Section 6. Some technical details and the time-dependent problem are presented in the 
Appendix. 

2. Linear Time-Independent Forecast/Analysis Cycle 

A general linear forecast/analysis cycle is 

y I = • 

yl = y{ + (y l - , 
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( 1 ) 

( 2 ) 



where the forecast, analysis and observation vectors at time step are denoted by y{, y£ and y° k 
respectively. The dynamics matrix advances the system from one analysis time to the next. 
The matrix H fc is the observation operator or generalized interpolation operator that maps forecast 
state variables to observation space. The dynamics and observation operator are taken to be linear 
and also time-independent for most of our results. This assumption on H* = H corresponds to 
a fixed, in situ, observing system. The assumption of linear dynamics (linearized about a non- 
linear trajectory) is not unreasonable for short times, i.e., up to one or two days in atmospheric 
models (Courtier and Talagrand, 1987). Stochastically forced linear models are also sometimes 
appropriate in an average sense (Blumenthal, 1991; Xue et al., 1994). The time-independence of 
the dynamics, = Ml, is a stronger restriction, although recent studies have shown that a sim- 
ple model linearized about the long-term mean flow and stochastically forced may be capable of 
reproducing aspects of observed storm-track variability (Whitaker and Sardeshmukh, 1998). In 
the Appendix we show how the time-independence requirements on the dynamics and observation 
operator can be relaxed. 

The state estimate y£ and the observations y° k are real vectors of length n and p respectively, 
usually with p <C n; 4** and H fc are real n x n and p x n matrices respectively; the gain K* is an 
n x p matrix. Typically in statistical interpolation and 3D-Var data assimilation methods the gain 
at time t/. is effectively given by 

K* = SiHj(H t S£Hj[ + Ri) _1 , (3) 

where S[ is a model of the forecast error covariance and R fc is a model of the observation error 
covariance (Cohn, 1997; Cohn et al., 1998). 

We assume the observation error to be additive and white in time with mean zero and p x p 
covariance matrix R*, and the forecast model error to be additive and white in time with mean zero 
and n x n covariance matrix Q™. In other words, denoting the true state of the atmosphere by y£, 
the observation vector y° k is related to the true state by 

yJ = H t y'+bJ, <bj) = 0 . (b£(b?f> = 5 kl R t , (4) 

and the evolution of the true state is given by 

yi = ^yt,+b?. <bD = o, (b?(bn T ) = W; < 5 > 
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b£ and b™ are the observation and forecast model error respectively; S k i is the Kronecker delta. 
Additionally the forecast model and observation errors are assumed to be uncorrelated with each 
other, (b£* (bf) T ) = 0, and with the initial error, (b™(y£ - y“) T ) = 0 and (b£(yo - y“) T ) = 0. 
We have neglected biases, which may be a particularly important part of the model error and should 
be estimated separately when present (Dee and da Silva, 1998). 

Under these assumptions, the time evolution of the analysis error y k = y t k — y£ is given by 

y k — Afcy^-i + b*; , (b/fe) = 0 , (b*bf) = SkiQk ) (yob J) = 0 , (6) 

where the system matrix f\ k and the system error covariance matrix Q k are defined by 

At = (I — Q t = (I - K t H t )Qr(l - K t H t ) r + K t R t K? ; (7) 

I is the n x n identity matrix. The system matrix and system error covariance matrix both depend 
on the gain matrix K*. The system matrix propagates analysis errors from one analysis time to 

the next. Roughly speaking, one expects analysis errors to grow in the forecast step (W fc ) and to 

be reduced by the analysis step (I - K fc H fc ). The stochastic forcing part b* of the analysis error 
propagation (6), with error covariance Q*, is due to both model and observation errors. 

Henceforth we assume time-independence of the dynamics and the observation operator, ty k = 
vp and H fc = H, and also of the model and observation error covariances, Q™ = Q m and = R. 
In this case, if the system (4), (5) is observable and if the (time-dependent) Kalman gain is used, 
then A fc — > A^ = A and Q fc — > Qoo — Q, and the eigenvalues of the steady-state system matrix 
A all lie inside the unit circle (Kalman, 1960). Suboptimal choices of the gain matrix may also 
yield an asymptotic steady state (e.g., Cohn and Todling, 1996). In any case, we assume that A^ 
and Qk tend to steady state, and that the steady-state system matrix A = Aoo is stable, that is, all 
its eigenvalues lie inside the unit circle. The steady-state analysis error covariance P is then given 
by the large-time limit (see the Appendix) 

P = lim (; y k yl ) , (8) 

and satisfies 

P = APA t + Q . (9) 

The linear matrix equation (9) is known as the discrete algebraic Lyapunov equation (DALE; Gajic 
and Qureshi, 1995). Some properties of its solution will be discussed in the next section. 
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3. Theory 

3.1 . Normal system matrix 


When the steady-state system matrix A is normal, the eigenvectors z* of A form a complete and 
orthonormal basis, and Q can be decomposed in this basis as 


n n 

Q = 12 ( z i Qz i) z * z i • ( 10 ) 

i=l jf=l 

The solution of (9) is then 


W it, it 

p = Ea‘Q(aT = EEt 

fc=0 1=1 7=1 i 


z*Qzj 


A i (A)A i (A) 




( 11 ) 


the notation A* (A) is used to denote the i-th eigenvalue of the matrix A, where the eigenvalues are 
ordered in decreasing magnitude; z* is the corresponding eigenvector of A and z* the conjugate 
transpose. Equation (11) shows that when A is normal, P is controlled by the proximity of the 
eigenvalues of A to the unit circle and by the relationship between the eigenvectors of A and Q. 
For instance, the lower bound 


a -< p >^ p *‘ = t^f (12) 

shows that if the leading eigenvalue of the steady-state system matrix, X x (A) , is near the unit circle, 
and if the projection z\Qz x of the steady-state system error covariance Q onto the corresponding 
eigenvector z x does not happen to be small, then P projects strongly onto the leading eigenvector 
Zx of the system matrix and the leading eigenvalue of P is large. In other words, the steady-state 
analysis error covariance P must have a large eigenvalue when the steady-state system matrix 
A has a weakly-damped eigenmode that is well-excited by the steady-state system error forcing. 
Alternatively, z^Qz x being large causes zJPzj to be large even when |Ax(A)| -C 1. 

The connection between the eigendecompositions of P and A is especially simple when the 
eigenvectors of A are also eigenvectors of Q (equivalently when A and Q commute), in which 
case Eq. (1 1) simplifies to 


p = £ 


Z*QZi 

1 — ! Ai (A) | 2 z * Zl ’ 


(13) 
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the eigenvalues of P are (z*Qzj)(l - |A;(A)| 2 ) -1 . In this case, when the leading eigenvalues of A 
are near the unit circle and the associated leading eigenvectors of A are well-excited by the system 
error, P will have large eigenvalues with eigenvectors that lie in the space spanned by the leading 
eigenvectors of A. 

In general, nonnormality of the dynamics and inhomogeneity of the observing network cause 
A to be nonnormal, and there is no simple connection between the eigenvectors of P and those of 
A. However, as will be shown in the next subsection, there is still a connection between the leading 
eigenvectors of P and more general stability properties of the steady-state system matrix A. 


3.2. Nonnormal system matrix 

When A is nonnormal its stability properties may not be well-described by its eigenvalues and 
eigenvectors (Trefethen et al., 1993; Farrell and Ioannou, 1996; Tippett et al., 2000). The damping 
of analysis errors may be much less efficient than predicted by the eigenvalues of A, and the 
dominant error structures may not be described well by eigenvectors of A. To obtain a more useful 
measure of the stability properties of A, we examine the quantity x r A iV y for any two real unit 
vectors x and y. This quantity gives the response in the x-direction to an analysis error in the 
y-direction after N forecast/analysis cycles in the steady-state regime. The maximum of |x T A iV y| 
defines the Euclidean matrix norm of A 77 , 



I T a AT 

max x A y 
IN| 2 =||y|| 2 =i 1 


(14) 


When A is normal and stable, jj A^ || 2 = ^A^zil = |Ai(A)| JV , and || A 7 ^ j| 2 decays nionotonically 
and exponentially as N increases since |Ai(A)| < 1. Therefore when A is normal and stable, 
HA^yfclk < llyfclh and the analysis error evolution equation (6) shows that near steady state, any 
transient analysis error growth must be due to the observational and model error forcing of the 
system error b^. In the general case when A is nonnormal, |x r A iV y | attains its largest value when 
x and y are respectively the leading eigenvectors of A' v ( A T ) N and (A T ) N A' v (equivalently, the 
leading left and right singular vectors of A^). The stability of A does not constrain [| A^ |j 2 to 
be monotonically decreasing as N increases. As a consequence, HA^y*^ need not be less than 
1 1 y A: 1 1 2 ; analysis errors in the steady-state regime can grow on transient time-scales even in the 
absence of observational and model error forcing. Asymptotically, however, jj A" v || 2 must behave 
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like |Ai(A)|‘ v (Horn and Johnson, 1985): 

lim exp -^7 log ||A iV || 2 = |A 1 (A)| . (15) 

(V— >oo iv 

and eventually analysis errors are reduced. In this limit of large N, the leading left and right singu- 
lar vectors of A v are approximately the leading eigenvectors of A and its adjoint A T respectively 
(Farrell, 1989). 

In the steady-state regime, the maximum response in the x-direction after N forecast/analysis 
cycles applied to a unit vector y is 

II^A^I^ = max |x T A"y| = J (x T A N ) (x T A Ar ) r . (16) 

l|y||2=i v 

A measure of the maximum response in the x-direction accumulated over all times-scales is ob- 
tained by summing the square of (16) over N : 


Y ||x t A"H! = ]Tx r A" (A T ) iV x = x r Bx , 


where the bound matrix B is defined by 


B = £a"(A’') w - 


Since B is a symmetric matrix, its eigenvalues have a “minimax” characterization (Golub and Van 


Loan, 1996), 


x t Bx 

Aj(B) = max min — - — . 

dim ( D)=i xG D X 1 X 


Therefore the eigenvectors of the bound matrix B order the directions in state space according to 
the maximum response in those directions. In this sense, the eigenvectors of B play a role similar to 
that of Lyapunov vectors, identifying directions in state space associated with accumulated growth. 

It is also useful to know which system error forcing produces the maximum or minimum 
steady-state analysis error variance since, as (7) shows, modifying the observing pattern can change 
the portion of the system error covariance due to model error. Applying the above arguments to 
the quantity || A^y j | 2 shows that the eigenvectors of the matrix B T given by (Tippett and March- 


esin, 1999a) 


b t = £( A i r A ‘ , 
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order the directions of state space according to the analysis error variance produced by system error 
forcing in those directions. When A is nonnormal, the eigenvectors of B and Bt are generally 
different, and therefore the directions of optimal response and forcing are different. 

3.3. Estimates derived from the bound matrix 

Given only limited information about the steady-state system error covariance Q, the bound ma- 
trix can be used to estimate the steady-state analysis error covariance P. Rigorous upper and 
lower matrix bounds for P are (Tippett and Marchesin, 1999a; see (54) of the Appendix for the 
corresponding bounds in the time-dependent case): 

A n (Q)B<P<A 1 (Q)B; (21) 

for two symmetric matrices X and Y, X < Y means that Y — X is positive semidefinite. When 
the ratio Ai(Q)/A n (Q) is near unity, these upper and lower bounds tightly restrict P. If more 
information about Q is available, tighter bounds can be obtained (see (56) of the Appendix and 
Tippett and Marchesin, 1999a). The estimates in (21) can be used to derive upper and lower 
bounds for the eigenvalues of P and therefore for the trace of P: 

A„(Q)Aj(B) < A*(P) < A 1 (Q)A i (B) , (22) 

and 

A n (Q)trB < trP < A x (Q)trB, / (23) 

where tr denotes the trace. Again, the upper and lower bounds are tight when the spectrum of Q 
is fairly flat. 

When the bounds in (21) are tight and P has a well-separated set of leading eigenvalues, the 
leading eigenvectors of P and B span approximately the same subspaces (Theorem 7.2.4 of Golub 
and Van Loan, 1996). However, if P does not have a well-separated set of leading eigenvalues, 
P — A!(Q)B may be small without the leading eigenvectors of P and B spanning approximately 
the same subspaces. 
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3.4. Estimates derived from approximation of the bound matrix 


For systems where n is large, calculating the bound matrix B directly may be as impractical as 
calculating P. However, if the series in (18) converges rapidly, a few terms may provide a good 
approximation for B, and therefore an additional theoretical tool. Iterative Lanczos methods can 
be used to calculate the leading eigenvalues and eigenvectors of the partial sum (Golub and Van 
Loan, 1996). The partial sum B L is defined by either 

Bz, = ]>>*(A T ) fc , (24) 

k - 0 

or recursively by B L = AB£_ X A T + I, B 0 = I. The cost of calculating eigenvectors of B L with an 
iterative method is about a factor L greater than that of iteratively calculating singular vectors of A. 
Monte Carlo simulations of (6) are another method of calculating B (see the Appendix). However, 
it is difficult to estimate a priori the number of terms needed to approximate B to a given accuracy. 
The relation 

l|B-B t || 2 =||A t + 1 B(A T ) I ' +1 || 2 <||A i+1 i||B|| 2 . (25) 

shows that the number of terms in (24) needed to approximate B to a given accuracy depends on 
how rapidly ||A L || 2 tends to zero. On the other hand, the size of ||B — B^|| 2 does not necessarily 
determine how well the eigenvectors of the truncated series B L approximate those of B, for even 
when || B — B L || 2 is not small, B and B^ may have approximately the same eigenvectors. For a 
normal system matrix, B and B L have precisely the same eigenvectors, independently of L > 1; 
the error patterns established after a single time step with the steady-state system matrix A are 
not subsequently altered. When A is nonnormal the situation is less clear, though there may be 
situations where the eigenvectors of B L approximate those of B for moderate values of L. We 
explore this point further in the numerical experiments of Section 5.2. 

The simplest approximation for the bound matrix (and one which allows some analytical treat- 
ment) is to take only the first two terms of the series in (18): 

BrjBj = I + AA t . (26) 

The eigenvalues of Bx are (1 + of (A)) and the eigenvectors are the left singular vectors of A; the 
singular values of A are defined by of ( A) = A;(AA r ). Since Bx < B, from (22) one has the lower 
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bound 


Ai(P) > AntQJAilBO = A„(Q)(1 + af(A)) , (27) 

and summing over i gives 

tr P > A n (Q)(l 4- tr AA t ) . (28) 

These two estimates highlight the role of nonmodal growth (ct* (A) > 1) on the size of the eigen- 
values of P. On the other hand, if A is normal, or nearly normal in the sense that nonmodal growth 
does not exceed the modal growth (tri(A) = |Ai(A)| < 1), bounds for the eigenvalues of P are 
(Mori et al., 1982) 

A - (Q) r^OT + Ai(Q) - Ai(p) - A ‘ (Q) + Al (q) t^Sp ' (29) 

When j Ax (A) | <C 1, as might be expected for an accurate and dense observing network for instance, 
(29) suggests that the eigenvalues of the steady-state analysis error covariance P are determined 
more by the steady-state system error covariance Q than by the steady-state system matrix A. 

4. Generalized advection model 

The analytical results presented above are illustrated in numerical experiments with two different 
dynamical models. We have chosen dynamics that are exponentially stable to demonstrate that 
even in the absence of modal instability, nonmodal growth is sufficient to cause the steady-state 
analysis error covariance to have dominant parts. The first dynamical model is a one-dimensional 
advection equation with a term added to cause nonmodal growth. Simple advection equations have 
been a valuable theoretical tool in data assimilation (e.g., Daley, 1992; Daley and Menard, 1993; 
Cohn, 1997; Mitchell and Daley, 1997). The simplicity of the dynamics allows the eigenvectors 
and singular vectors to be computed analytically. 

4.1. Model 

The generalized one-dimensional advection equation considered is 

(j, t + a/j, x + c' (x) / j. = 0, (30) 
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with initial and periodic boundary conditions 


ix(x,t = 0) = Ho(x), /j.{x = 0,£) = n(x = 1, i) , (31) 

respectively. The quantity a plays the role of a constant zonal velocity and —c'(x) that of a zonally 
varying vertical velocity. When the vertical velocity is identically zero, c'(x) = 0, the dynamics is 
normal. As we shall show, the larger the vertical velocity, the more nonnormal is the dynamics. 

The physical basis for the model is the following. Considering an xz — plane representing the 
zonal and vertical directions, a two-dimensional nondivergent wind field can be written in terms 
of a stream function ip(x,z) as — i/j z x + ip x z. Suppose that the stream function has the form 
ip(x, z) = —az — c(x)\ a is the constant zonal velocity; the vertical velocity is —d{x) where ' is 
used to indicate the derivative with respect to x. The streamlines of such a wind field are shown in 
Figure 1 for c(x) taken to be 

c(x) = ^|(1 + cos 2tt(z — 0.5)) 4 , (32) 

with c 0 = 0.15 and a = 1/9. These non-dimensional parameters correspond to a zonal length scale 
of 3105 km, a vertical length scale of 1 km, a zonal velocity of 15 m/s, a maximum vertical velocity 
of 10 cm/s and a time-scale of 6.4 hours. Suppose that a quantity p is advected by such a wind 
field and that initially the 2 -dependence of p has the form p(x, z, t = 0) = ix{x, t = 0) exp(— z). 
Then the equation for the time evolution of /i is (30). 

The operator U* T that advances the solution r time units is given by 

t)} = n(x — ar,t) exp [(c(x — ar) — c(x))/a ] . (33) 

An “exact” discrete dynamics is obtained by restricting U/ T to act on functions defined on a spa- 
tially regular grid with spacing Ax such that the Courant number arj Ax is an integer. If the time- 
step t is small compared to the advection time-scale 1/a, the factor exp[(c(z — ar) — c(x))/a] 
can be approximated by exp (— d(x)r) and in a single time-step the solution is seen to grow in 
regions where the vertical velocity is positive ( —c'(x ) > 0). This behavior is consistent with the 
physical interpretation of p having vertical density stratification that decreases with height. This 
growth is not evident from the eigenvalues of W r . The eigenvalues of U^ T are exp (— 2-Karli) with 
corresponding eigenfunctions exp [—c{x)/a + 2 nlxi], for any integer l. Regardless of the vertical 
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velocity —c'(x), all the eigenvalues of 4f r have unit modulus. Only in the case that c'(x ) = 0 
and hence U* r is normal, does this property actually imply that solutions neither decay nor grow in 
time. 


4.2. Singular values and vectors 

The effect of applying the operator Uf T once is explained by its singular values (or more generally 
s-numbers; Gohberg and Krein, 1969). The singular values of 4> r are the square roots of the 
eigenvalues of v|/ T \|/* (or of 4>*4J X ). The adjoint of equation (30) with respect to the usual inner 
product on the space of square integrable functions is 

ix t - a/ijj + c'(x)fjL = 0 . (34) 

Hence, the adjoint operator 4** is 

ty*{jj,(x, £)} = fi(x + ar, t ) exp [— (c(x + ar) — c(:r))/a] . (35) 

The eigenvalues of 4** are exp (2tv arli) with corresponding eigenfunctions exp [c{x)/a 4- 2nlxi\. 
Using the definitions in (33) and (35), 

tU T M/*{^(a:)} = fx{x ) exp [2(c(a; — ar) — c(x))/a ] . (36) 

The operator is a multiplication operator and its spectrum consists of the values taken on by 

the multiplying factor exp [2(c(r — ar) — c(x))/a\) for 0 < x < 1 (Halmos, 1967, Problem 52).* 

Although the eigenvalues of the operator 4* r are necessarily the complex conjugates of those of 

its adjoint 4/*, the eigenfunctions of 4* r and 4>* can be very different, leading to transient growth 

and limiting predictability (Farrell, 1989). If c / ( x ) = 0, then 4/ T is normal, and the eigenfunctions 

of Ur r and 4/* are the same. More generally the cosine S of the angle between the eigenfunctions 

*In general, the spectrum of an operator can be divided into three disjoint sets: the point spectrum, the residu- 
al spectrum and the continuous spectrum. For a multiplication operator A induced by a multiplier 4>: (i) the point 
spectrum is equal to the set of complex numbers for which has positive measure; (ii) the residual spectrum is 

empty; and (iii) the continuous spectrum is the essential range of (p minus the point spectrum (Halmos, 1967, Problem 
66). An eigenvalue A is in the point spectrum of A if and only if there exists a function / such that Af = A/. For an 
eigenvalue A in the continuous spectrum there is no eigenfunction but there is a sequence of approximate eigenfunc- 
tions fj such that \\Afj — A/j || — > 0. Therefore, for non-constant c(x) the operator UJ T U** has a continuous spectrum 
and no eigenfunctions but does have approximate eigenfunctions, for example those given by the eigenvectors of the 
discrete approximation 
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of vl/ r and 4** is given by 


^2 = f o J Q exp (2 (c(x) - c(y))/a) dxdy . (37) 

When the vertical velocity — c'(x) is large, the cosine S will be small, indicating that the eigen- 
functions of 4/ r and 4/* are nearly orthogonal. Similarly, the smaller the zonal velocity a is, the 
smaller S is. The quantity S is also a measure of the sensitivity of the eigenvalues to perturbations 
of the operator 4* r (Golub and Van Loan, 1996, Section 7.2.2). 

The singular values of the spatially discrete dynamics 4^ are obtained by restricting the values 
of x to be grid points. That is, according to Eq. (36), the n singular values a m of 4^ are (not 
ordered by size): 

= exp [(c(a: m - ar) - c(a; m ))/a] , x m = (m- l)/n, m — l, ...,n. (38) 

This expression shows that if ar <C 1, the size of the largest singular value is an (exponentially) 
increasing function of the maximum vertical velocity and the time-step r. In Fig. 2, c(x) and 
s(x) = exp [(c(x — ar) — c(x))/a] are plotted for c(x) given by (32) and r = 0.5, 1.0, 2.0. For 
larger values of r the magnitude of the largest singular value increases and its localization moves 
upstream. 

The n x n matrices and 4^*4^ are diagonal, a reflection of the fact that 4> r 4>* and 

4** 4> T are multiplication operators; the squared singular values cr^ appear along the diagonal. The 
diagonal entries of 4^*4*^ are those of 4^4^* shifted by am units, the Courant number, as can 
be seen by comparing the expression for 4**4* r with that of given in (36). Thus the left 

singular vector of corresponding to a m is an n-vector with value 1 at the m-th position and 
zero elsewhere, and the right singular vector corresponding to the same singular value has 1 at the 
(m — am) position and is elsewhere zero. When there are repeated singular values, the associated 
singular vectors of are not uniquely defined; only the subspace they span is uniquely defined. 
The localized structure of the singular vectors seen here is observed, to some extent, in more 
physically realistic models (e.g., Buizza et al., 1997). The left and right singular vectors of 4^ are 
different in general but for this problem it is only their ordering that is different. 
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4.3. Kalman filter experiments 


Three sets of experiments were performed. In all experiments, the solution of the DALE is obtained 
by iterating the time-dependent Kalman filter to steady state and using the resulting steady-state 
gain K to compute the system matrix. In the first set of experiments we take d{x) = 0 and 
the dynamics is normal; an inhomogeneous observing pattern is the source of system matrix 
nonnormality (Daley, 1992). In the second and third sets of experiments nonnormal dynamics 
is considered and c(x ) is that given in (32). The second and third sets differ in their choices of 
observing pattern. In the second set of experiments, the nine trailing left singular vectors of the 
dynamics are observed, leading to a system matrix with large singular values. In the third set of 
experiments, the nine leading left singular vectors of the dynamics are observed. This is equivalent 
to observing upstream (experiment 2) and in the middle (experiment 3) of the instability (Todling 
and Ghil, 1996). Observations are assimilated at unit time intervals, r = 1; the r subscript and the 
d superscript on 4^ will be dropped from the notation. A grid spacing of Ax = 1/36 is used, so 
that n = 36. In all cases, the observations are made at grid points and the number of observations 
p — 9 is one fourth of the total number of grid points. 

The observation error is taken to be spatially homogeneous and uncorrelated, R = 0.025I P for 
all experiments; l p is the p x p identity matrix. We consider two types of model error covariances: 
one with Q m = 0.0251, for which errors at different grid points are uncorrelated, and another 
whose entries are constructed as 

Q 7 j =* ij + P(d ij ;L), ' (39) 

where dij is the distance between the i-th and j- th grid points and p is the compactly supported 
piecewise-rational correlation function with correlation length L = 0.2 (Gaspari and Cohn, 1999, 
Eq. 4.10). The trace of Q m is normalized to 0.025n. A 36 x 36 “checkerboard plot” of the elements 
of the Q m given by Eq. (39) is shown in Fig. 3. 

In the first experiment c'(x) = 0 and the dynamics is normal. Observations are taken at the 
first 9 grid points. Figure 4a shows the singular values and the absolute values of the eigenvalues 
of the steady-state system matrix A for Q m = 0.0251. Though 4* is normal, the inhomogeneous 
observing pattern makes A nonnormal; the difference between the singular values and magnitudes 
of the eigenvalues is one indication of nonnormality. The singular values of the system matrix 


17 



control how analysis errors are modified by one forecast/analysis cycle in the steady-state regime. 
Figure 4a shows that there is a nine-dimensional subspace where analysis errors are reduced after a 
single forecast/analysis cycle; the magnitude of analysis errors outside of this subspace remain un- 
changed. The eigenvalues of A are all less than one in modulus, demonstrating that asymptotically 
all analysis errors are reduced. 

The spectra of the observation and model error covariance matrices R = 0.025 l p and Q m = 
0.0251 are flat by construction. However, the spectrum of P shown in Fig. 4b varies by an order 
of magnitude between its largest and smallest eigenvalues (note the logarithmic scale on panel 
b). This feature is similar to that seen in the example of Tippett and Marchesin (1999 b) where 
for the same dynamics, assimilation of a single observation produced a steady-state analysis error 
covariance whose largest and smallest eigenvalues differed by a factor equal to the ratio of the 
domain size to the advection speed; here, the ratio of the domain size to the advection speed is 9. 
The clustering of the eigenvalues of P in groups of four is the result of the spatial homogeneity 
of the model and observational errors and of the advection speed, 4 grid-points/unit time. Also 
shown in Fig. 4b are the upper and lower estimates for the eigenvalues of P obtained from (22). 
The corresponding results (not shown) for the experiment with spatially correlated model error are 
qualitatively similar to the spatially uncorrelated model error case though the bounds from (22) are 
not tight because A n (Q) <C A X (Q) for the spatially correlated case. 

The magnitude of the eigenvalues and singular values of *P for the nonnormal case are shown 
in Fig. 5. The singular values of U> are the values of the function plotted (dashed line) in Fig. 2, 
evaluated at the grid points. As was shown analytically, though the eigenvalues of the dynamics 
all have magnitude one, the presence of singular values greater than unity allows transient growth 
to occur. To examine the effect of not observing the growing modes of the dynamics, first only 
the trailing 9 left singular vectors of *P were observed. We take H = [u n _ p+ i, . . . , u„] T , where 
u m is the m-th left singular vector of 'P and p = 9. Figure 2 shows that this choice of H corre- 
sponds to observing to the left or upstream of the instability. Figure 6a shows that as in the case 
with normal dynamics, the analysis errors are reduced in a nine-dimensional subspace by a single 
forecast/analysis cycle. Since the observations were chosen to lie in the region where errors are 
not amplified by the dynamics, the smallest singular values of A shown in Fig. 6a are smaller than 
those in the normal case, shown in Fig. 4a. Not observing the growing modes of the dynamics 
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leads to A, like W, having 13 singular values greater than unity, indicating that there is a thirteen- 
dimensional subspace associated with transiently growing analysis errors. The asymptotic rate at 
which analysis errors are reduced, as given by the eigenvalues of A (crosses in Fig. 6a), is roughly 
the same as in the normal case (crosses in Fig. 4a). Figure 6b shows that not observing grow- 
ing modes of the dynamics leads to P having eigenvalues that are larger than in the normal case 
(compare Fig. 4b). Similar results (not shown) were found for the spatially correlated model error 
case. 

Figure 7a shows that with the same nonnormal dynamics, if the 9 leading left singular vectors 
of the dynamics are observed (H = [ux, . . . , u p ] T ), the singular values of A are reduced and most 
of the growing modes of A are removed. This H corresponds to observing directly in the middle 
of the instability. The result of this change in observing pattern is a reduction of the analysis error 
to a level comparable to that of the normal case (compare Fig. 7b and Fig. 4b). This reduction in 
size of the eigenvalues of P is suggested by the inequality in (27). However, as the relation in (27) 
is a lower bound on the eigenvalues of P, reducing the singular values of the steady-state system 
matrix does not necessarily reduce the asymptotic analysis error levels. If the analysis error levels 
were strictly controlled by the singular values of the system matrix, then the optimal observing 
strategy would be simply to minimize the singular values of the system matrix. We return to this 
point in Section 5 and in the Conclusions. 

We have shown how the eigenvalues of the bound matrix B are related to those of the steady- 
state analysis error covariance P. Now we examine how the eigenvectors of B project onto those 
of P. For the three pairs of experiments described, the invariant subspaces (subspaces spanned by 
the eigenvectors) of P and B were compared in the following manner: let W = [wi, . . . , w n ] and 
Z = [zx, . . . , z n ] be the n x n matrices whose columns are the eigenvectors, in order of decreasing 
eigenvalue, of P and B respectively. Consider the matrix M defined to be the product M = Z T W. 
The entries of M are M i: , = zjw j. That is, the i-th row of M contains the components of the 
z-th eigenvector of B written in the basis given by the eigenvectors of P. Gelaro et al. (1998) used 
a similar approach, referring to M as a projection matrix, to compare the subspaces spanned by 
initial and final singular vectors. If P and B have the same eigenvectors with the same ordering, 
then M = I; if only the ordering is different then M is the identity matrix with some columns 
permuted. Since W and Z are orthogonal matrices, the entries of M have absolute value between 
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zero and one. 

The structure of M for the three cases considered is shown in Figs. 8 - 10. In these “checker- 
board” plots of the elements of M, the darkness of the ij - th square is proportional to the absolute 
value of M, r The darkest squares indicate vectors that are parallel; white squares indicate vectors 
that are orthogonal. The image on the left is for Q m = 0.0251 and that on the right for the spatially 
correlated Q m . Figs. 8-10 (left panels) show that in all the cases when the model error covariance 
is proportional to the identity matrix, P and B have the same eigenvectors; only the ordering is dif- 
ferent. Additionally in all the cases the difference in ordering is not very great. When the spatially 
correlated Q m is used, the eigenvectors of P and B do not coincide (Figs. 8-10, right panels). 
However, some diagonal dominance is still seen in the plots of M, indicating that the leading in- 
variant subspace of B is still a good approximation of that of P; that is, the leading eigenvectors of 
P project mostly onto the leading eigenvectors of B. 

5. Damped baroclinic wave model 

The next illustration uses the dynamics of a model of damped baroclinic waves (Farrell, 1989). 
Similar baroclinic models have been studied extensively in the context of optimal growth and 
predictability (Farrell, 1985; Farrell and Ioannou, 1993). 

5.1. Model 

Using the same scaling parameters as in Farrell (1989), the non-dimensionalized potential vorticity 
equation based on the non-dimensional stream function 

$(x, y, z, t ) = ip(z, t)e z / 2 e lkx cos (ly) (40) 

(41) 

(42) 


is 


^- + ikz 
dt 


tfzz -[ of +-] ll> 


+ ik(fl + 1 )ib = 0 . 


with boundary conditions 


9 (ipz + t ) ~ **(! ” W = 0 , 2 = 0, 


dt 
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(i’z + ^0 - zA:^ = 0, 2 = z T , (43) 

where a = sjk 2 + l 2 and T is the non-dimensionalized Ekman damping parameter taken here to 
be 0.2z, corresponding to a verticaleddy viscosity of 4.5m 2 /s; (3 = 0.533 and a unit time interval is 
equivalent to 9.3 hours. The k = l =. 2 wave with meridional wavelength 3100 km is considered. 
The upper boundary condition at z T = 4.0 corresponds to a height of 40 km. 

Discretizing equation (41) and the boundary conditions (42) and (43) using centered differences 
on a regular grid in the vertical coordinate 2 gives a system of equations of the form 

^ F jmlpm , (44) 

m~l 

where Vj(t) = -fi{zj,t), Zj is the j- th level, and n is the number of levels; details of F jm can be 
found in Farrell (1989). The matrix that advances the solution by time r is U^ T = exp (rF). 
Figure 1 1 shows the singular values and eigenvalues of W T with the number of levels n = 36, for 
r = 0.6452 corresponding to a time step of 6 hours, which is used here. The qualitative features of 
the singular values and eigenvalues of the baroclinic dynamics are similar to those of the general- 
ized advection equation shown in Fig. 5, with eigenvalues that have approximately unit magnitude 
and several singular values greater than unity, indicating transient, nonmodal growth. The pertur- 
bation that leads to maximum initial growth was calculated in Farrell (1989) and has most of its 
structure concentrated near the surface. In the absence of Ekman damping, the baroclinic wave 
dynamics has one exponentially increasing and one exponentially decaying eigenmode in addition 
to the neutral eigenmodes. The vertical dissipation used here is sufficient to damp the exponential- 
ly growing mode (Farrell, 1989). However, an effect of the relatively low vertical resolution used 
here is that the exponential growth rate is quite small but nonzero. 

5.2. Kalman filter experiments 

The steady-state system matrix A and analysis error covariance P are calculated as in Section 4.3. 
The observation error is taken to be uncorrelated, with R = 0.1lj>; the forecast model error given in 
(39) is used with L — 5 km and the normalization tr Q m /n = 1. To give an indication of the effect 
of changing the observing pattern, the squared Frobenius norm of A fc , j|A fc |||. .= tr A fc (A t )*\ is 
plotted as a function of k in Fig. 12 for three different observing patterns. This quantity can be 
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interpreted as the expected value (y^y k) of the response of the unforced system y fc+l = Ay fc to 
a spatially uncorrelated random initial condition y 0 with (yoYo ) = l‘» that is, (yjyk) — ||A fc ||^, 
and initially ||A°||^ = ||l|||* = n = 36. This unforced, time-independent system, like the bound 
matrix, can be used as a tool to study the steady state of the time-dependent system (6). The 
three observing patterns are (1) observations of ip at every fourth vertical level; (2) observations 
of the leading 6 right singular vectors of the dynamics; and (3) observations of the leading 6 left 
singular vectors of the dynamics. Both the leading left and the leading right singular vectors of 
the dynamics have most of their structure localized near the surface. Figure 12 shows that with 
observing pattern 1, 1| [||r has a large initial transient in comparison to that of observing patterns 

2 and 3. This behavior is a reminder that the asymptotic gain K, while optimal in the asymptotic 
steady state, can give far from optimal performance in the unsteady regime. After the initial peak, 
||A fe |||p decays faster for pattern 1 than for patterns 2 and 3. The size of the steady-state analysis 
error variance depends on the behavior of the steady-state system matrix on all time-scales. In 
particular, from (18), 

oo 

trB = ^||A fc |||, " (45) 

k = o 

and therefore the trace of the bound matrix is proportional to the area under the graph of ||A*|||. 
plotted as a function of k. Although the temporal behavior of the steady-state system matrices 
associated with the three observing patterns is different, the areas under the graphs in Fig. 12 are 
not so different, and thus the steady-state values of the total analysis error variance are comparable; 
for pattern 1, tr P = 236.4, for pattern 2, tr P = 251.7, and for pattern 3, tr P = 309.2. 

The remainder of the results shown use observing pattern 1. Figure 13a shows two properties of 
the steady-state system matrix A that lead to large asymptotic analysis error levels: (i) the system 
matrix has eigenvalues inside but near the unit circle; and (ii) it has many singular values greater 
than unity. Note that the largest singular values of the system matrix are larger than those of the 
dynamics, another indication of the suboptimality of the asymptotic gain outside of the steady- 
state regime; we show later that these growing modes of the system matrix are “unlikely” in the 
steady state since they project mostly onto the trailing eigenvectors of the steady-state analysis 
error covariance. Figure 13b shows that the steady-state analysis error covariance has large eigen- 
values with the bounds from (22) predicting well the shape but not the values of the analysis error 
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covariance spectrum. 

We examine the eigenvectors of the steady-state analysis error covariance in the same manner 
used in Section 4.3. Figure 14 shows checkerboard plots of the elements of M = Z* W for W taken 
to be the eigenvectors of P, and for Z taken to be (a) the eigenvectors of B, (b) the eigenvectors 
of Bi (left singular vectors of A) and (c) the right singular vectors of A. Figure 14(a) shows that 
the leading eigenvectors of the bound matrix B are nearly the same as those of the asymptotic 
analysis error covariance matrix. Figure 14(b) shows that the leading invariant subspace of Bi 
(leading left singular vectors of A) is also a good approximation for that of the asymptotic analysis 
error covariance matrix. The dark squares in the lower left-hand comer of Fig. 14(c) show that the 
trailing right singular vectors of the steady-state system matrix are a good approximation for the 
leading invariant subspace of the asymptotic analysis error covariance matrix. Additionally, the 
leading right singular vectors of the system matrix project strongly onto the trailing eigenvectors 
of P (see upper right-hand comer of Fig. 14(c)). The leading right and left singular vectors of the 
steady-state system matrix A are therefore approximately orthogonal since they approximate the 
trailing and leading eigenvectors, respectively, of P, and since the leading and trailing eigenvectors 
of P are orthogonal by construction. 

Similar relations between right and left singular vectors of dynamics have been noted else- 
where. For the discrete dynamics of Sec. 4, we saw that the left and right singular vector corre- 
sponding to any given singular value are orthogonal, being different columns of the identity matrix. 
Gelaro et al. (1998) observed near-orthogonality of leading right and left singular vectors of the 
dynamics in a more realistic model. An explanation for this similar behavior of the steady-state 
system matrix is that since the system matrix is stable, its leading left singular vectors must project 
mostly onto its trailing right singular vectors, which are orthogonal to its leading right singular 
vectors by construction, so that growth is not maintained. 

We stated earlier that the eigenvectors of B, by ordering state space according to accumulated 
growth, play a role analogous to that of Lyapunov vectors. That analogy is supported by the 
relations observed here among the eigenvectors of B, the singular vectors of the steady-state system 
matrix, and the eigenvectors of the steady-state analysis error covariance matrix. Szunyogh et al. 
(1997) noted that generally leading right singular vectors of the dynamics are quite different from 
leading Lyapunov vectors, while leading left singular vectors resemble leading Lyapunov vectors. 
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Also, Swanson et al. (1998), observing all state variables, and using perfect-model, 4D-Var data 
assimilation, found that the dominant analysis errors were well described by the Lyapunov vectors. 

6. Summary and Conclusions 

The properties of a general linear, time-independent forecast/analysis system were investigated. 
The steady-state analysis error covariance matrix of such a system satisfies a discrete algebraic 
Lyapunov equation (DALE). Using properties of the DALE we examined why the steady-state 
analysis error covariance might have a dominant part and how that dominant part can be approxi- 
mated. Nonnormality of the system matrix is a key feature; the system matrix is the linear operator 
that combines the forecast and analysis steps. If the system matrix is normal, the properties of the 
DALE and hence of the steady-state analysis error covariance are trivially related to the eigenval- 
ues and eigenvectors of the system matrix. The role of observation and forecast model errors is 
also important since their presence prevents repeated iterations of the forecast/analysis cycle from 
relaxing to modal behavior. However, the view taken here has been to assume that the impact of the 
system matrix on the steady-state analysis error covariance is greater than that of the observation 
and forecast model error forcing. This choice is a practical one since the forecast model error is 
poorly known. 

The bound matrix provides a quantitative connection between the structure of the steady-state 
analysis error covariance and the nonmodal stability properties of the system matrix. The eigen- 
values of the bound matrix give estimates for the eigenvalues of the steady-state analysis error 
covariance matrix, and in some circumstances the leading eigenvectors of the bound matrix ap- 
proximate those of the steady-state analysis error covariance matrix. This bound matrix estimate 
approach extends naturally to the time-dependent problem. The eigenvectors of the bound ma- 
trix, identifying the directions in state space associated with accumulated past growth, play a role 
similar to that of Lyapunov vectors in deterministic systems. Methods of computing the leading 
singular values and vectors of the system matrix might be adapted to compute the leading eigen- 
vectors of an approximate bound matrix. In the simplest approximation, the leading eigenvectors 
of the bound matrix are approximated by the leading left singular vectors of the system matrix, 
and the leading eigenvalues of the bound matrix depend on the squares of the singular values of 
the system matrix. The singular values and vectors of the system matrix are in general different 
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from those of the dynamics. If growing modes of the dynamics, in particular leading left singular 
vectors of the dynamics, are not observed, then the analysis step will not be able to reduce the 
associated errors in the first guess. In that case, growing modes of the system matrix are related to 
those of the dynamics. 

Failure to observe the leading left singular vectors of the dynamics was seen to lead to a system 
matrix with large singular values and elevated steady-state analysis error variances in a numerical 
example using a generalized advection equation. This result has implications for targeted observ- 
ing strategies. However, there are several caveats. The first is the simplicity of the dynamics, 
in particular, its time-independence. Second, in this work the focus has been on the analysis er- 
ror covariance, whereas investigations of observing strategies as well as initial ensemble selection 
typically have used forecast skill as the measure of optimality (Gelaro et al., 1998). The choice of 
whether to observe the modes that grew (left singular vectors) or the modes that will grow (right 
singular vectors) depends on whether forecast or analysis accuracy is the optimality criterion. If 
analysis accuracy is desired, then observations of the leading left singular vectors of the dynamics 
are generally more important than observations of the leading right singular vectors of the dynam- 
ics. The reason for this conclusion is that if observations of the leading right singular vectors of 
the dynamics are used, the size of the forecast error that projects onto the leading left singular 
vectors of the dynamics depends on the observation error amplified by the forecast dynamics. The 
nonnormality of the dynamics implies that the subspaces spanned by the leading right and left sin- 
gular vectors may be nearly orthogonal. In that case, with observations of only the leading right 
singular vectors, the analysis is unable to reduce the first guess error that projects onto the leading 
left singular vectors of the dynamics. On the other hand, if observations of the left singular vectors 
are used, the first guess error projecting onto the leading left singular vectors of the dynamics is 
reduced by the observations, and the size of the resulting analysis error depends on the observa- 
tion error unamplified by the forecast dynamics. Of course, applying the same arguments to the 
steady-state forecast error covariance leads to the conclusion that observations of the leading right 
singular vectors of the dynamics reduce singular values of the forecast error propagation system 
matrix and consequently forecast error. 

Since the dependence of the steady-state analysis error covariance on the singular values of the 
system matrix presented in this work is in the form of a lower bound, increasing the singular values 
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of the system matrix guarantees that the steady-state analysis error variance will be increased, but 
reducing them does not guarantee reduction of the variance. That is to say, observing leading the 
left singular vectors of the dynamics may not be the optimal strategy since the steady-state analysis 
error covariance does not depend on a single application of the system matrix. More generally, the 
dependence of the bound matrix on the system matrix should be used to examine the role of the 
observing pattern in determining steady-state analysis error levels. It must be remembered that the 
observing pattern also plays a role in the system error covariance matrix, which includes terms for 
both the model and the observation error covariances. Failure to observe the leading eigenvectors 
of the model error covariance leads to large system error covariance eigenvalues and consequently 
to high analysis error levels. In a similar fashion, the dependence of the bound matrix on the 
dynamics gives a means of investigating covariance evolution schemes based on low-dimensional 
representation of the dynamics. 

In summary we find that: 

• Nonnormality of the forecast/analysis system is a key factor in determining the properties of 
the steady-state analysis error covariance. 

• Theoretical estimates based on the bound matrix show the dependence of the steady-state 
analysis error covariance on the dynamics and observing pattern through the system matrix. 

• In the simplest approximation, the dominant eigenvectors of the bound matrix are approxi- 
mated by the left singular vectors of the system matrix; when growing modes of the dynamics 
are not observed, they are related to the left singular vectors of the dynamics. 

This last point seems to have implications in the context of approximating error covariances 
via ensemble methods (e.g., Evensen, 1994; Houtekamer et al., 1996), at least for linear systems. 
The accuracy of an ensemble method depends on capturing growing modes of the system and on 
representing the analysis uncertainty well. The work here suggests that it may not always be possi- 
ble to satisfy these criteria simultaneously. Choosing perturbations that correspond to leading right 
singular vectors, while capturing the growing modes, would not necessarily result in an ensemble 
whose sample covariance approximates well the analysis error covariance. In fact, in the model 
for damped baroclinic waves, the leading eigenvectors of the steady-state analysis error covariance 
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projected strongly onto the trailing right singular vectors of the system matrix. Conversely, the 
sample forecast error covariance resulting from a small initial ensemble that approximates well the 
analysis uncertainty may not explain well forecast errors due to growing modes. Possible solutions 
are to project analysis errors on to the growing modes or to calculate the growing modes weighted 
by the analysis error (Barkmeijer et al., 1998; Swanson et al., 1998). 
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Appendix 

Define P* = (y kYk) where y k is defined by (6), and suppose that (6) is in steady state: A fc = A 
and Q k = Q. Then 

Pyk = ((Ay/;_ i + bjfc)(Ayfc_! + b*;) T ) 

, = A (yk-iYk-i) A T + Q + A (y fc -ibj[) + ( b k yl ’_ x ) A T . 

But (yjfc-xb^) = 0 since 

h — 1 

y/c-i = A t - 1 yo + ^ Ak " i "’ b y. < 47 > 

j=l 

(y 0 bfc) = 0, and (b^b^) = 0 fory ^ k. Therefore, P k = AP fc _!A T -I- Q, or 

p fc = A*P 0 (A T ) fc + J2 AjQ ( AT )' ' (48) 

i - o 

Since the eigenvalues of A lie inside the unit circle, the first term in (48) tends to zero, and the 
terms of the series in (48) decay geometrically and the series converges uniformly (Golub and Van 
Loan, 1996, Lemma 7.3.2). Therefore, the limit 

P = lim Pfc (49) 

k—too 


27 



exists and is the solution of P = APA T + Q. 

The bound matrix approach can be extended to the time-dependent case as follows. Again 


define P k = (y k yj) where y k is defined by (6). Then using the arguments above, 

Pfc = AfcPjfc-! A^ + Qfc , (50) 

and we can define a time-dependent bound matrix B k by the iteration 

Bfc = AfcBfc-iA^ + 1 , (51) 

with B 0 = I. Subtracting (50) from (51) multiplied by a scalar j k gives 

(TfcBjt — Pfc) = Afc( 7 fcB*_i — P/;_i)Aj + ( 7*1 — Qfc) . (52) 

Define the nonnegative constants 7 ^ and j k by 

7 fc = max(Ai(P 0 ), max Ai(Qj)) , and 7 ^ = min(A„(P 0 ), min A n (Qj)) . (53) 


Since ( 7 * I — Q fc ) is negative semidefinite for all k, it follows from (52) that ( 7 ^ B fc — P fc ) is negative 
semi definite for all k. Likewise it follows that (j k B k — Pfc) is positive semidefinite for all k. In 
other words, 

7 * B fc < P* < 7 ^B fc , (54) 

with the relation in ( 21 ) following as a special case. 

If there is some information available about the system error covariance, it can be used in the 
definition of the bound matrix. Suppose we have a time-dependent model Q k for the system error 
covariance Q k and time-independent constants 7 “ and 7 + such that 7 _ Qfc < Q k < 7 + Qfc. Then 
defining a bound matrix B^ by 

Bjt = AfcBfc-tA^ + Qfc , (55) 

it follows that 

7 “Bfc < Pfc < 7 + Bfc. (56) 

If the system error covariance Qfc is known at least approximately, Monte Carlo methods can be 
used to calculate Pfc by approximating the expectation with an ensemble average (Evensen, 1994). 



Such an approach is similar to that of system simulation (Houtekamer et al., 1996) where an attempt 
is made to account for all sources of error; it is also conceptually similar to the “breeding of 
growing modes” (BGM) method (Toth and Kalnay, 1993). In the BGM method the factor (I — 
KfcHfc) appearing in the system matrix is approximated by a constant or by a “regional rescaling” 
factor, and there is no explicit stochastic forcing. 

In the steady-state problem, the dependence of the properties of B on the nonnormality of 
the system matrix highlights the role played by nonnormality in stochastically forced systems. If 
the forcing term is not included, the iteration in (6) with A& = A is the power method without 
renormalization (Golub and Van Loan, 1996, Section 7.3.1). For this reason, it has been stated 
that for linear, stationary dynamics the BGM method is similar to the power method and that in 
this case, the leading BGM vector is the leading eigenvector of the dynamics (Buizza and Palmer, 
1995). If the system matrix A is normal, this conclusion is correct even in the presence of stochastic 
forcing, since in that case, the leading eigenvector of B is identical to the leading eigenvector of 
A. However, when the system matrix is nonnormal, the eigenvectors of the bound matrix B are in 
general not eigenvectors of A. Another reason not to interpret BGM as the power method is that 
the power method may be excessively slow to converge for nonnormal matrices (Golub and Van 
Loan, 1996, Lemma 7.3.1). 

The eigenvectors of the bound matrix B do not correspond to Lyapunov vectors of the unforced 
system. For time-independent dynamics the method for computing Lyapunov vectors (Legras 
and Vautard, 1996, Section 5) reduces to orthogonal iteration, a standard method for calculating 
the subspace spanned by the leading eigenvectors (Golub and Van Loan, 1996, Section 7.3.2); 
the Lyapunov vectors themselves are not eigenvectors but an orthogonal basis for the invariant 
subspaces of A (for an alternative Lyapunov vector definition see Trevisan and Pancotti, 1998). 
However, conceptually the eigenvectors of B do play the role of Lyapunov vectors of the system 
with spatially uncorrelated, homogeneous random forcing, since they order the directions of state 
space according to accumulated growth. 
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Fig. 1. The streamlines of ip{x, z) = —az — c(x) for c(x) given by (32). 


35 




X 


Fig. 2. Solid line c(x); s(:r) = exp [(c(x — ar) — c(x))/a] for r = 0.5 (dotted line), r 
(dashed line), and r = 2.0 (dotted-dashed line). 
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Fig. 3. Checkerboard plot of the elements of the spatially correlated Q m used in the experiments 
and grey scale. 
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(b) 




Fig. 4. Singular values and magnitude of the eigenvalues of the steady-state system matrix for 
normal dynamics ( c'(x ) = 0) and observations at the first 9 grid points. Panel (a) shows the 
singular values (plus signs) and magnitude of the eigenvalues (crosses) of A for Q m = 0.0251. 
Panel (b) shows the eigenvalues of P (crosses), and upper and lower bounds from .(22) (dotted 
lines). 
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Fig. 5. Singular values (plus signs) and magnitude of the eigenvalues (crosses) of 4* for the non- 
normal dynamics case. 
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(a) (b) 




Fig. 6. As in Fig. 4, but for nonnormal dynamics and an observation operator that corresponds to 
observing the trailing left singular vectors of the dynamics M*. 
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Fig. 8. Checkerboard plots representing the projection of the eigenvectors of P onto the eigenvec- 
tors of B for the normal dynamics case. The image on the left is for Q m = 0.0251 and that on the 
right for the spatially correlated model error. 
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Fig. 9. As in Fig. 8, but for nonnormal dynamics and with observations of the trailing left singular 
vectors of the dynamics. 
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Fig. 10. As in Fig. 9, but with observations of the leading left singular vectors of the dynamics. 
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Fig. 11. Singular values (plus signs) and magnitude of the eigenvalues (crosses) of for the 
damped baroclinic wave model. 
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Fig. 12. The quantity ||A fc |||. plotted as a function of k for the damped baroclinic wave dynamics 
and observing patterns pattern 1 (diamonds), pattern 2 (plus signs) and pattern 3 (asterisks) as 
described in the text. 
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Fig. 13. Singular values and magnitude of the eigenvalues of the system matrix for dynamics of 
the damped baroclinic wave model. Panel (a) shows the singular values (plus signs) and magnitude 
of the eigenvalues (crosses) of A. Panel (b) shows the eigenvalues of P (crosses), and upper and 
lower bounds from (22) (dotted lines). 
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Fig. 14. Checkerboard plots of the elements of the matrix M = Z*W for Z taken to be (a) the 
eigenvectors of B, (b) the eigenvectors of B L and (c) the right singular vectors of A. 
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